Correlated geminal wave function for molecules: an efficient resonating valence bond 

approach 

Michele CasulaQ Claudio AttaccaliteQ and Sandra Sorelkd 
International School for Advanced Studies (SISSA) Via Beirut 2,4 34014 Trieste , 
Italy and INFM Democritos National Simulation Center, Trieste, Italy 
(Dated: February 2, 2008) 

We show that a simple correlated wave function, obtained by applying a Jastrow correlation 
term to an Antisymmetrized Geminal Power (AGP), based upon singlet pairs between electrons, is 
particularly suited for describing the electronic structure of molecules, yielding a large amount of the 
correlation energy. The remarkable feature of this approach is that, in principle, several Resonating 
Valence Bonds (RVB) can be dealt simultaneously with a single determinant, at a computational cost 
growing with the number of electrons similarly to more conventional methods, such as Hartree-Fock 
(HF) or Density Functional Theory (DFT). Moreover we describe an extension of the Stochastic 
Reconfiguration (SR) method, that was recently introduced for the energy minimization of simple 
atomic wave functions. Within this extension the atomic positions can be considered as further 
variational parameters, that can be optimized together with the remaining ones. The method is 
applied to several molecules from Li2 to benzene by obtaining total energies, bond lengths and 
binding energies comparable with much more demanding multi configuration schemes. 



I. INTRODUCTION 

The comprehension of the nature of the chemical bond 
deeply lies on quantum mechanics; since the seminal work 
by Heitler and London very large steps have been 
made towards the possibility to predict the quantitative 
properties of the chemical compounds from a theoretical 
point of view. Mean field theories, such as HF have been 
successfully applied to a wide variety of interesting sys- 
tems, although they fail in describing those in which the 
correlation is crucial to characterize correctly the chemi- 
cal bonds. For instance the molecular hydrogen H2, the 
simplest and first studied molecule, is poorly described by 
a single Slater determinant in the large distance regime, 
which is the paradigm of a strongly correlated bond; in- 
deed, in order to avoid expensive energy contributions - 
the so called ionic terms - that arise from two electrons of 
opposite spin surrounding the same hydrogen atom, one 
needs at least two Slater determinants to deal with a spin 
singlet wave function containing bonding and antibond- 
ing molecular orbitals. Moreover at the bond distance it 
turns out that the resonance between those two orbitals 
is important to yield accurate bond length and binding 
energy, as the correct rate between the ionic and covalent 
character is recovered. Another route that leads to the 
same result is to deal with an AGP wave function, which 
includes the correlation in the geminal expansion; Barbi- 
ellini in Ref.0gave an illuminating example of the beauty 
of this approach solving merely the simple problem of the 
H 2 molecule. 

On the other hand the variational methods based on 
the Configuration Interaction (CI) technique, which is 
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able to take into account many Slater determinants, have 
been shown to be successful for small molecules (e.g. Be? 
Q)- In these cases it is indeed feasible to enlarge the 
variational basis up to the saturation, the electron cor- 
relation properties are well described and consequently 
all the chemical properties can be predicted with ac- 
curacy. However, for interesting systems with a large 
number of atoms this approach is impossible with a rea- 
sonable computational time. Coming back to the H? 
paradigm, it is straightforward to show that a gas with 
N H-2 molecules, in the dilute limit, can be dealt accu- 
rately only with 2 N Slater determinants, otherwise one 
is missing important correlations due to the antibonding 
molecular orbital contributions, referred to each of the N 
H2 molecules. Therefore, if the accuracy in the total en- 
ergy per atom is kept fixed, a Cl-likc approach does not 
scale polynomially with the number of atoms. Although 
the polynomial cost of these Quantum Chemistry algo- 
rithms - ranging from N 5 to N 7 - is not prohibitive, a 
loss of accuracy, decreasing exponentially with the num- 
ber of atoms is always implied, at least in their simplest 
variational formulations. This is related to the loss of size 
consistency of a truncated CI expansion. On the other 
hand, this problem can be overcome by coupled cluster 
methods, that however in their practical realization are 
not variational^. 

An alternative approach, not limited to small 
molecules, is based on DFT. This theory is in principle 
exact, but its practical implementation requires an ap- 
proximation for the exchange and correlation function- 
als based on first principles, like the Local Density ap- 
proximation (LDA) and its further gradient corrections 
(GGA), or on semi empirical approaches, like BLYP and 
B3LYP. For this reason, even though much effort has 
been made so far to go beyond the standard functionals, 
DFT is not completely reliable in those cases in which 
the correlation plays a crucial role. Indeed it fails in de- 
scribing HTc superconductors and Mott insulators, and 
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in predicting some transition metal compounds proper- 
ties, whenever the outermost atomic d-shell is near-half- 
filled, as for instance in the high potential iron proteins 
(HiPIP)||. Also i?2 molecule in the large distance regime 
must be included in that list, since the large distance 
Born-Oppenheimer energy surface, depending on Van der 
Waals forces, is not well reproduced by the standard func- 
tionals, although recently some progress has been made 
to include these important contributions 0. 

Quantum Monte Carlo (QMC) methods are alternative 
to the previous ones and until now they have been mainly 
used in two versions: 

• Variational Monte Carlo (VMC) applied to a wave 
function with a Jastrow factor that fulfills the cusp 
conditions and optimizes the convergence of the CI 
basis H El; 

• Diffusion Monte Carlo (DMC) algorithm used to 
improve, often by a large amount, the correlation 
energy of any given variational guess in an auto- 
matic manner [9j. 

Hereafter we want to show that a large amount of 
the correlation energy can be obtained with a single de- 
terminant, using a size-consistent AGP- Jastrow (J AGP) 
wave function. Clearly our method is approximate and in 
some cases not yet satisfactory, but in a large number of 
interesting molecules we obtain results comparable and 
even better than multi determinants schemes based on 
few Slater determinants per atom that are affordable by 
QMC only for rather small molecules. 

Moreover, we have extended the standard SR method 
to treat the atomic positions as further variational pa- 
rameters. This improvement, together with the possi- 
bility to work with a single determinant, has allowed 
us to perform a structural optimization in a non triv- 
ial molecule like the benzene radical cation, reaching the 
chemical accuracy with an all-electron and feasible vari- 
ational approach. 

The paper is organized as follows: In Sec. [H] we in- 
troduce the variational wave function, that is expanded 
over a set of non orthogonal atomic orbitals both in the 
determinantal AGP and the Jastrow part. This basis set 
is consistently optimized using the method described in 
Sec. IHII that. as mentioned before, allows also the geom- 
etry optimization. Results and discussions are presented 
in the remaining sections. 



II. FUNCTIONAL FORM OF THE WAVE 
FUNCTION 

In this paper we are going to extend the application 
of the JAGP wave function, already used to study some 
atomic systems ^(J. We generalize its functional form 
in order to describe the electronic structure of a generic 
cluster containing several nuclei. With the aim to deter- 
mine a variational wave function, suitable for a complex 



electronic system, it is important to satisfy, as we re- 
quire in the forthcoming chapters, the "size consistency" 
property: if we smoothly increase the distance between 
two regions A and B each containing a given number 
of atoms, the many-electron wave function factorizes 
into the product of space-disjoint terms ^ = ^a^I^b 
as long as the interaction between the electrons coupling 
the different regions A and B can be neglected. In this 
limit the total energy of the wave function approaches 
the sum of the energies corresponding to the two space- 
disjoint regions. This property, that is obviously valid 
for the exact many-electron ground state, is not always 
fulfilled by a generic variational wave function. 

Our variational wave function is defined by the prod- 
uct of two terms, namely a Jastrow J and an antisym- 
metric part = J^agp)- If t nc former is an explicit 
contribution to the dynamic electronic correlation, the 
latter is able to treat the non dynamic one arising from 
near degenerate orbitals through the geminal expansion. 
Therefore our wave function is highly correlated and it is 
expected to give accurate results especially for molecular 
systems. The Jastrow term is further split into a two 
body and a three body factors, J = J2J3, described in 
the following subsections after the AGP part. 



A. Pairing determinant 

As is well known, a simple Slater determinant provides 
the exact exchange electron interaction but neglects the 
electronic correlation, which is by definition the missing 
energy contribution. In the past different strategies were 
proposed to go beyond Hartree-Fock theory. In particu- 
lar a sizable amount of the correlation energy is obtained 
by applying to a Slater determinant a so called Jastrow 
term, that explicitly takes into account the pairwise in- 
teraction between electrons. QMC allows to deal with 
this term in an efficient wav|ll|. On the other hand, 
within the Quantum Chemistry community AGP is a well 
known ansatz to improve the HF theory, because it im- 
plicitly includes most of the double-excitations of an HF 
state. 

Recently we proposed a new trial function for atoms, 
that includes both the terms. Only the interplay between 
them yields in some cases, like Be or Mg, an extremely 
accurate description of the correlation energy. In this 
work we extend this promising approach to a number of 
small molecular systems with known experimental prop- 
erties, that are commonly used for testing new numerical 
techniques. 

The major advantage of this approach is the inclusion 
of many CI expansion terms with the computational cost 
of a single determinant, 

that allow us to extend the calculation with a full struc- 
tural optimization up to benzene, without a particularly 
heavy computational effort on a single processor machine. 
For an unpolarizcd system containing N electrons (the 
first N/2 coordinates are referred to the up spin elec- 
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trons) the AGP wave function is a x pairing matrix 
determinant, which reads: 

*AGp(ri,...,rjv) = det ($ AGP (r t ,r J+N/2 )) , (1) 

and the gcminal function is expanded over an atomic ba- 
sis: 

<^ GP (rV)= ]T A^> a , ; (rT)^. m (r^), (2) 

l 7 m,a,b 

where indices m span different orbitals centered on 
atoms a, b, and are coordinates of spin up and down 
electrons respectively. The geminal functions may be 
viewed as an extension of the simple HF wavefunction, 
based on molecular orbitals, and in fact the geminal func- 
tion coincide with HF only when the number M of non 
zero eigenvalues of the A matrix is equal to N/2. Indeed 
the general function [2] can be written in diagonal form 
after an appropriate transformation: 

M 

<S> AGP (r\rt) = J2^Mr V )Mr l ), (3) 

k 

where <pk(j) = Ylj a L JL k.j,a4>j.a{'£) are just the molecular 
orbitals of the HF theory whenever M = N/2. Notice 
that with respect to our previous pairing function for- 
mulation also off-diagonal elements are now included in 
the A matrix, which must be symmetric in order to de- 
fine a singlet spin orbital state. Moreover it allows one to 
easily fulfill other system symmetries by setting the ap- 
propriate equalities among different A;, m . For instance 
in homo- nuclear diatomic molecules, the invariance un- 
der reflection in the plane perpendicular to the molecular 
axis yields the following relation: 

= (_l)Pm+PnX b > a (4) 

where p m is the parity under reflection of the m— th or- 
bital. 

Another important property of this formalism is the 
possibility to describe resonating bonds present in many 
structures, like benzene. A A^ b „ different from zero rep- 
resents a chemical bond formed by the linear combination 
of the m-th and n-th orbitals belonging to a-th and 6-th 
nuclei. It turns out that resonating bonds can be well de- 
scribed through the geminal expansion switching on the 
appropriate A°; b „ coefficients: the relative weight of each 
bond is related to the amplitude of its A. 

Also the spin polarized molecules can bc treated within 
this framework, by using the spin generalized version of 
the AGP (GAGP), in which the unpaired orbitals are ex- 
panded as well as the paired ones over the same atomic 
basis employed in the geminal |l2j|. As already mentioned 
in the introduction of this chapter, the size consistency 
is an appealing feature of the AGP term. Strictly speak- 
ing, the AGP wave function is certainly size consistent 
when both the compound and the separated fragments 
have the minimum possible total spin, because the gem- 
inal expansion contains both bonding and antibonding 



contributions, that can mutually cancel the ionic term 
arising in the asymptotically separate regime. Moreover 
the size consistency of the AGP, as well as the one of the 
Hartree-Fock state, holds in all cases in which the spin 
of the compound is the sum of the spin of the fragments. 
However, similarly to other approaches 3, for spin polar- 
ized systems the size consistency does not generally hold, 
and, in such cases, it may be important go beyond a single 
AGP wave function. Nevertheless we have experienced 
that a single reference AGP state is able to describe ac- 
curately the electronic structure of the compound around 
the Born-Oppenheimer minimum even in the mentioned 
polarized cases, such as in the oxygen dimer. 

The last part of this section is devoted to the nuclear 
cusp condition implementation. A straightforward cal- 
culation shows that the AGP wave function fulfills the 
cusp conditions around the nucleus a if the following lin- 
ear system is satisfied: 

(Is. 2s) 

E KifoA* = ' = ^K;io,^:r = R„i. 

3 c,3 

(5) 

for all b and j'; in the LHS the caret denotes the spherical 
average of the orbital gradient. The system can be solved 
iteratively during the optimization processes, but if we 
impose that the orbitals satisfy the single atomic cusp 
conditions, it reduces to: 

£ A J c ;f<MRa) = 0, (6) 

c(#a),j 

and because of the exponential orbital damping, if the 
nuclei are not close together each term in the previous 
equations is very small, of the order of exp(— |R a — R c |). 
Therefore, with the aim of making the optimization 
faster, we have chosen to use Is and 2s orbitals satisfy- 
ing the atomic cusp conditions and to disregard the sum 
© in Eq. [SJ In this way, once the energy minimum is 
reached, also the molecular cusp conditions JSJ arc rather 
well satisfied. 



B. Two body Jastrow term 

As it is well known the Jastrow term plays a crucial 
role in treating many body correlation effects. One of the 
most important correlation contributions arises from the 
electron-electron interaction. Therefore it is worth us- 
ing at least a two-body Jastrow factor in the trial wave 
function. Indeed this term reduces the electron coales- 
cence probability, and so decreases the average value of 
the repulsive interaction. The two-body Jastrow function 
reads: 

( N \ 

J 2 (ri,...,rjv) = exp ^!i(r y ) , (7) 
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where u(rij) depends only on the relative distance = 
— Tjl between two electrons and allows to fulfill the 
cusp conditions for opposite spin electrons as long as 
u i r ij) ~ * ~2~ f° r smau electron-electron distance. The 
pair correlation function u can be parametrized success- 
fully by few variational parameters. We have adopted 
two main functional forms. The first is similar to the one 
given by Ceperley [l3[ : 

u(r) = ^(l-e- r / F ), (8) 

with F being a free variational parameter. This form 
for u is particularly convenient whenever atoms are very 
far apart at distances much larger than F, as it allows 
to obtain good size consistent energies, approximately 
equal to the sum of the atomic contributions, without 
changing the other parts of the wave function with an 
expensive optimization. Within the functional form JSJl, 
it is assumed that the long range part of the Jastrow, 
decaying as a power of the distance between atoms, is 
included in the 3-body Jastrow term described in the 
next subsection. The second form of the pair function 
m, particularly convenient at the chemical bond distance, 
where we performed most of the calculations, is the one 
used by Fahy Q: 



with a different variational parameter b. 

In both functional forms the cusp condition for antipar- 
allcl spin electrons is satisfied, whereas the one for paral- 
lel spins is neglected in order to avoid the spin contamina- 
tion. This allows to remove the singularities of the local 
energy due to the collision of two opposite spin electrons, 
yielding a smaller variance and a more efficient VMC cal- 
culation. Moreover, due to the Jastrow correlation, an 
exact property of the ground state wave function is re- 
covered without using many Slater determinants, thus 
considerably simplifying the variational parametrization 
of a correlated wave function. 



C. Three Body Jastrow term 

In order to describe well the correlation between elec- 
trons the simple Jastrow factor is not sufficient. Indeed 
it takes into account only the electron-electron separa- 
tion and not the individual electronic position and r^. 
It is expected that close to nuclei the electron correla- 
tion is not accurately described by a translationally in- 
variant Jastrow, as shown by different authors, see for 
instance Ref. For this reason we introduce a fac- 

tor, often called three body (electron-electron-nuclcus) 
Jastrow, that explicitly depends on both the electronic 
positions and Tj . The three body Jastrow is chosen to 
satisfy the following requirements: 

• The cusp conditions set up by the two-body Jas- 
trow term and by the AGP are preserved. 



• It does not distinguish the electronic spins other- 
wise causing spin contamination. 

• Whenever the atomic distances are large it factor- 
izes into a product of independent contributions lo- 
cated near each atom, an important requirement to 
satisfy the size consistency of the variational wave 
function. 

Analogously to the pairing trial function in Eq. [21 we 
define a three body factor as: 

J 3 (n,...,r N ) = exp ^^^(r^rjoj 

Mri,r 3 -) = E 9ti^aA^bM, (io) 

l,m,a,b 

where indices I and m indicate different orbitals located 
around the atoms a and b respectively. Each Jastrow 
orbital ip a< i(r) is centered on the corresponding atomic 
position R a . We have used Gaussian and exponential or- 
bitals multiplied by appropriate polynomials of the elec- 
tronic coordinates, related to different spherical harmon- 
ics with given angular momentum, as in the usual Slater 
basis. Analogously to the geminal function $agp, when- 
ever the one particle basis set {tp a ,i} is complete the ex- 
pansion ltTU|) is also complete for the generic two particle 
function $j(r,r'). In the latter case, however, the one 
particle orbitals have to behave smoothly close to the 
corresponding nuclei, namely as: 

Vv( r ) - Vw(Ra) ^ |r-R a | 2 , (11) 

or with larger power, in order to preserve the nuclear 
cusp conditions 

For the s-wavc orbitals we have found energetically 
convenient to add a finite constant cj/(JV— 1). As shown 
in the Appendix [5] a non zero value of the constant q 
for such orbitals ip a ,i is equivalent to include in the wave 
function a size consistent one body term. As pointed out 
in Ref. Il6l it is easier to optimize a one body term im- 
plicitly present in the 3-body Jastrow factor, rather than 
including more orbitals in the determinantal basis set. 

The chosen form for the 3-body Jastrow (|10f> is simi- 
lar to one used by Prendergast et al. and has very 
appealing features: it easily allows including the symme- 
tries of the system by imposing them on the matrix gf'^ 
exactly as it is possible for the pairing part (e.g. by re- 
placing with g%i b n in Eq. |3}. It is size consistent, 
namely the atomic limit can be smoothly recovered with 
the same trial function when the matrix terms gf'^ for 
a =/= b approach zero in this limit. Notice that a small non 
zero value of g?' for a ^ b acting on p-wave orbitals can 
correctly describe a weak interaction between electrons 
such as the the Van dcr Waals forces. 



5 



III. OPTIMIZATION METHOD 



Sk 4 is easily expressed in terms of these operators: 



We have used the Stochastic Reconfiguration (SR) 
method already described in Ref . 0, that allows to min- 
imize the energy expectation value of a variational wave 
function containing many variational parameters in an 
arbitrary functional form. The basic ingredient for the 
stochastic minimization of the wave function deter- 
mined by p variational parameters is the 
solution of the linear system: 



p 

fe=0 



s j>k Aa k = (*\O k (AI - H)\V) 



(12) 



where the operators O k are defined on each N electron 
configuration x = {ri, . . . , r^} as the logarithmic deriva- 
tives with respect to the parameters a k : 



o k (x) 



d 
da k 



ln^fx) for k > 0, 



(13) 



while for k = Ok is the identity operator equal to one 
on all the configurations. The (p + 1) x (p + 1) matrix 



J 



(14) 



and is calculated at each iteration through a standard 
variational Monte Carlo sampling ; the single iteration 
constitutes a small simulation that will be referred in the 
following as "bin" . After each bin the wave function pa- 
rameters are iteratively updated (ctk — > oik + Ao^/Acno), 
and the method is convergent to an energy minimum for 
large enough A. Of course for particularly simple func- 
tional form of ^(x), containing e.g. only linear CI co- 
efficients, much more efficients optimization schemes do 
exist El. 



SR is similar to a standard steepest descent (SD) 
calculation, where the expectation value of the energy 
E(a k ) = ^I^P is optimized by iteratively changing the 
parameters a, according to the corresponding derivatives 
of the energy (generalized forces): 



fk 



dE 
da k 



(■*\O k H + HO k + (d ak H)\V) , (*|0 fc |*)(¥|ff|*) 



(15) 



namely: 



form similar to the steepest descent: 



OL k 



a k + At/ fc . 



(16) 



At is a suitable small time step, which can be taken fixed 
or determined at each iteration by minimizing the energy 
expectation value. Indeed the variation of the total en- 
ergy AE at each step is easily shown to be negative for 
small enough At because, in this limit 



AE = -At ^ f- + 0(At 2 ). 



Thus the method certainly converges at the minimum 
when all the forces vanish. Notice that in the defini- 
tion of the generalized forces (|15|) we have generally as- 
sumed that the variational parameters may appear also in 
the Hamiltonian. This is particularly important for the 
structural optimization since the atomic positions that 
minimize the energy enter both in the wave function and 
in the potential. 

In the following we will show that similar consider- 
ations hold for the SR method, that can be therefore 
extended to the optimization of the geometry. Indeed, 
by eliminating the equation with index k = from the 
linear system (I12I1 . the SR iteration can be written in a 



k 



where the reduced p x p matrix s is: 

Sj,k = s j,k — Sj,0S0,fc 
and the At value is given by: 

At 



2(A 



<*!*> 



(17) 



(18) 



(19) 



^From the latter equation the value of At changes during 
the simulation and remains small for large enough energy 
shift A. However, using the analogy with the steepest 
descent, convergence to the energy minimum is reached 
also when the value of At is sufficiently small and is kept 
constant for each iteration. Indeed the energy variation 
for a small change of the parameters is: 



AE = -AtJ2s;Jfifj 



It is easily verified that the above term is always negative 
because the reduced matrix s, as well as s _1 , is positive 
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definite, being s an overlap matrix with all positive eigen- 
values. 

For a stable iterative method, such as the SR or the 
SD one, a basic ingredient is that at each iteration the 
new parameters a' are close to the previous a accord- 
ing to a prescribed distance. The fundamental difference 
between the SR minimization and the standard steepest 
descent is just related to the definition of this distance. 
For the SD it is the usual one defined by the Cartesian 
metric A Q = J^k \ a 'k ~ a k\ 2 , instead the SR works cor- 
rectly in the physical Hilbert space metric of the wave 
function yielding A a = ^,j( a i ~ a i)i a 'j — <Xj), 
namely the square distance between the two normalized 
wave functions corresponding to the two different sets of 
variational parameters {a'} and {afc}. Therefore, from 
the knowledge of the generalized forces fk , the most con- 
venient change of the variational parameters minimizes 
the functional Ai? + AA Q , where AE is the linear change 
in the energy AE = — J^i fi{ a i~ a i) an d A is a Lagrange 
multiplier that allows a stable minimization with small 
change A Q of the wave function "J. The final iteration 
HI 711 is then easily obtained. 

The advantage of SR compared with SD is obvious 
because sometimes a small change of the variational pa- 
rameters correspond to a large change of the wave func- 
tion, and the SR takes into account this effect through 
the Eq. 1171 In particular the method is useful when a 
non orthogonal basis set is used as we have done in this 
work. Indeed by using the reduced matrix s it is also 
possible to remove from the calculation those parame- 
ters that imply some redundancy in the variational space. 
As shown in the Appendix^ a more efficient change in 
the wave function can be obtained by updating only the 
variational parameters that remain independent within 
a prescribed tolerance, and therefore, by removing the 
parameters that linearly depend on the others. A more 
stable minimization is obtained without spoiling the ac- 
curacy of the calculation. A weak tolerance criterium 
e ~ 10 -3 , provides a very stable algorithm even when the 
dimension of the variational space is large. For a small 
atomic basis set, by an appropriate choice of the Jastrow 
and Slater orbitals, the reduced matrix s is always very 
well conditioned even for the largest system studied, and 
the above stabilization technique is not required. Instead 
the described method is particularly important for the 
extension of QMC to complex systems with large num- 
ber of atoms and/or higher level of accuracy, because in 
this case it is very difficult to select - e.g. by trial and 
error - the relevant variational parameters, that allow a 
well conditioned matrix s for a stable inversion in i|17|) . 



A. Setting the SR parameters 

In this work, instead of setting the constant A, we 
have equivalently chosen to determine At by verifying 
the stability and the convergence of the algorithm at 
fixed At value, which can be easily understood as an 



inverse energy scale. The simulation is stable when- 
ever I /At > A cut , where A cllt is an energy cutoff that 
is strongly dependent on the chosen wave function and 
is generally weakly dependent on the bin length. When- 
ever the wave function is too much detailed, namely has 
a lot of variational freedom, especially for the high en- 
ergy components of the core electrons, the value of A cllt 
becomes exceedingly large and too many iterations are 
required for obtaining a converged variational wave func- 
tion. In fact a rough estimate of the corresponding num- 
ber of iterations P is given by PAt » 1/G, where G is 
the typical energy gap of the system, of the order of few 
electron Volts in small atoms and molecules. Within the 
SR method it is therefore extremely important to work 
with a bin length rather small, so that many iterations 
can be performed without much effort. 

In a Monte Carlo optimization framework the forces fk 
are always determined with some statistical noise r]k, and 
by iterating the procedure several times with a fixed bin 
length the variational parameters will fluctuate around 
their mean values. These statistical fluctuations are sim- 
ilar to the thermal noise of a standard Langevin equation: 

d t a k = fk + Vk, (20) 

where 

(Vk(t)Vk'{t')) = 2T noise 6(t - i>M'- ( 21 ) 

The variational parameters oik, averaged over the 
Langevin simulation time (as for instance in Fig^ for 
t > 2H~ 1 ), will be close to the true energy minimum, but 
the corresponding forces fk = —d ak E will be affected by 
a bias that scales to zero with the thermal noise T no i se , 
due to the presence of non quadratic terms in the energy 
landscape. 

Within a QMC scheme, one needs to estimate T no i se , 
by increasing the bin length as clearly T no i se oc 
1/Bin length, because the statistical fluctuations of the 
forces, obviously decreasing by increasing the bin length, 
are related to the thermal noise by Eq. Thus there is 
an optimal value for the bin length, which guarantees a 
fast convergence and avoid the forces to be biased within 
the statistical accuracy of the sampling. 

An example is shown in Fig. ^ for the optimization of 
the Be atom, using a DZ basis both for the geminal and 
the three-body Jastrow part. The convergence is reached 
in about 1000 iteration with At = 0.005i? _1 . However, 
in this case it is possible to use a small bin length, yield- 
ing a statistical accuracy in the energy much poorer than 
the final accuracy of about 0.05mH. This is obtained by 
averaging the variational parameters in the last 1000 iter- 
ations, when they fluctuate around a mean value, allow- 
ing a very accurate determination of the energy minimum 
which satisfies the Euler conditions, namely with fk = 
for all parameters. Those conditions have been tested by 
an independent Monte Carlo simulation about 600 times 
longer than the bin used during the minimization. As 
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FIG. 1: Example of the convergence of the SR method for the variational parameters of the Be atom, as a function of the 
number of stochastic iterations. In the upper(lower) panel the Jastrow (geminal) parameters are shown. For each iteration, a 
variational Monte Carlo calculation is employed with a bin containing 15000 samples of the energy, yielding at the equilibrium a 
standard deviation of ~ 0.001877. For the first 200 iteration At = 0.00125JT" 1 , for the further 200 iterations At = 0.0025H -1 , 
whereas for the remaining ones At — 0.005i/ _1 . 



shown in Fig. [21 the Euler conditions are fulfilled within 
statistical accuracy even when the bin used for the mini- 
mization is much smaller than the overall simulation. On 
the other hand if the bin used is too small, as we have 
already pointed out, the averaging of the parameters is 
affected by a sizable bias. 
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Whenever it is possible to use a relatively small bin in 
the minimization, the apparently large number of iter- 
ations required for equilibration does not really matter, 
because a comparable amount of time has to be spent 
in the averaging of the variational parameters, as shown 
in Fig. H The comparison shown in the Ref. ^3 about 
the number of the iterations required, though is clearly 
relevant for a deterministic method, is certainly incom- 
plete for a statistical method, because in the latter case 
an iteration can be performed in principle in a very short 
time, namely with a rather small bin. 

It is indeed possible that for high enough accuracy the 
number of iterations needed for the equilibration becomes 
negligible from the computational point of view. In fact 
in order to reduce, e.g. by a factor ten, the accuracy in 
the variational parameters, a bin ten times larger is re- 
quired for decreasing the thermal noise T noise by the same 
factor, whereas a statistical average 100 times longer is 
indeed necessary to reduce the statistical errors of the 
variational parameters by the same ratio. This means 
that the fraction of time spent for equilibration becomes 
ten times smaller compared with the less accurate simu- 
lation. 



FIG. 2: Calculation of the derivative of the energy with re- 
spect to the second Z in the 2p orbital of the geminal function 
for the Be atom. The calculation of the force was obtained, 
at fixed variational parameters, by averaging over 10 7 sam- 
ples, allowing e.g. a statistical accuracy in the total energy 
of 0.07m//. The variational parameters have been obtained 
by an SR minimization with fixed bin length shown in the 
x label. The parameter considered has the largest deviation 
from the Euler conditions. 



B. Structural optimization 

In the last few years remarkable progresses have been 
made to develop Quantum Monte Carlo (QMC) tech- 
niques which are able in principle to perform struc- 
tural o ptim ization of molecules and complex systems 
|20t l2lL |22|. Within the Born-Oppcnhcimcr approxima- 
tion the nuclear positions can be considered as further 
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variational parameters included in the set {cti} used for 
the SR minimization l|17|) of the energy expectation value. 
For clarity, in order to distinguish the conventional varia- 
tional parameters from the ionic positions, in this section 
we indicate with {ci} the former ones, and with R; the 
latter ones. It is understood that R" = a&, where a 
particular index k of the whole set of parameters {at} 
corresponds to a given spatial component (y = 1, 2, 3) of 
the i— th ion. Analogously the forces (|15|l acting on the 
ionic positions will be indicated by capital letters with 
the same index notations. 

The purpose of the present section is to compute 
the forces F acting on each of the M nuclear positions 
{Ri, . . . , Rm}, being M the total number of nuclei in the 
system: 



F(R a ) = -VR a £({cJ,R a ), 



(22) 



with a reasonable statistical accuracy, so that the itera- 
tion (|17J) can be effective for the structural optimization. 
In this work we have used a finite difference operator 
for the evaluation of the force acting on a given 
nuclear position a: 



F(R„) 



E 



E(Kg + ARq) - E(Rg - ARq) 

2AR 



(23) 



0(AR 2 



where AR a is a 3 dimensional vector. Its length AR 
is chosen to be 0.01 atomic units, a value that is small 
enough for negligible finite difference errors. In order to 
evaluate the energy differences in Eq. 1231 we have used 
the space- warp coordinate transformation [23l |2^ | briefly 
summarized in the following paragraphs. According to 
this transformation also the electronic coordinates r will 



be translated in order to mimic the right displacement of 
the charge around the nucleus a: 



r, = r, 



AR„ LU a (Ti), 



where 



w a (r) 



F(|r-R„|) 



(24) 



(25) 



F(r) is a function which must decay rapidly; here we 
used F(r) = pc as suggested in Ref. l24l 

The expectation value of the energy depends on AR, 
because both the Hamiltonian and the wave function de- 
pend on the nuclear positions. Now let us apply the 
space- warp transformation to the integral involved in the 
calculation; the expectation value reads: 



B(R+AR) = /^(^ r ^-»- r 

/drJ AR (r)^ R (r(r)) 



, R (r(r))£f R (r(r)) 



(26) 



where J is the Jacobian of the transformation and here 
and henceforth we avoid for simplicity to use the atomic 
subindcx a. The importance of the space warp in re- 
ducing the variance of the force is easily understood for 
the case of an isolated atom a. Here the force acting on 
the atom is obviously zero, but only after the space warp 
transformation with u> a = 1 the integrand of expression 
(|26[) will be independent of AR, providing an estimator 
of the force with zero variance. 

Starting from Eq. |5HJ it is straightforward to explicitly 
derive a finite difference differential expression for the 
force estimator, which is related to the gradient of the 
previous quantity with respect to AR, in the limit of the 
displacement tending to zero: 



J 



F(R) = -( lim 
v ' \arh AR 



El) + 2((H)( ^ torfjV**)) _ {H |A Hm q ]og(jV»*))) , 



|AR|- 



(27) 



r 



where the brackets indicate a Monte Carlo like average 

over the square modulus of the trial wave function, 

is the finite difference derivative as defined in (|23|l . and 



E, 



{%\H\* 

WW) 



is the local energy on a configuration x 
where all electron positions and spins are given. In anal- 
ogy with the general expression (|15f) of the forces, we 
can identify the operators corresponding to the space- 
warp change of the variational wave function: 



0k = a^ 1os(j ^* ar) 



(28) 



The above operators l|28|l are used also in the definition 
of the reduced matrix s for those elements depending on 
the variation with respect to a nuclear coordinate. In 



this way it is possible to optimize both the wave function 
and the ionic positions at the same time, in close analogy 
with the Car-Parrinello[25| metho d ap plied to the min- 
imization problem. Also Tanaka |26j tried to perform 
Car-Parrincllo like simulations via QMC, within the less 
efficient steepest descent framework. 

An important source of systematic errors is the depen- 
dence of the variational parameters Ci on the ionic config- 
uration R, because for the final equilibrium geometry all 
the forces fi corresponding to Cj have to be zero, in or- 
der to guarantee that the true minimum of the potential 
energy surface (PES) is reached |2^. As shown clearly 
in the previous subsection, within a QMC approach it is 
possible to control this condition by increasing systcmat- 
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ically the bin length, when the thermal bias T no i se van- 
ishes. In Fig. |21 we report the equilibrium distance of the 
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FIG. 3: Plot of the equilibrium distance of the Lii molecule 
as a function of the inverse bin length. The total energy 
and the binding energy are reported in Tables [HI and II I II 
respectively. The triangles (full dots) refer to a simulation 
performed using 1000 (3000) iterations with At = 0.015//~1 
(At = 0.005//" 1) and averaging over the last 750 (2250) itera- 
tions. For all simulations the initial wavefunction is optimized 
at Li — Li distance 6 a.u. 

Li molecule as a function of the inverse bin length, so that 
an accurate evaluation of such an important quantity is 
possible even when the number of variational parameters 
is rather large, by extrapolating the value to an infinite 
bin length. However, as it is seen in the picture, though 
the inclusion of the 3s orbital in the atomic AGP basis 
substantially improves the equilibrium distance and the 
total energy by ~ ImH, this larger basis makes our simu- 
lation less efficient, as the time step At has to be reduced 
by a factor three. 

We have not attempted to extend the geometry op- 
timization to the more accurate DMC, since there are 
technical difficulties , and it is computationally much 
more demanding. 

C. Stochastic versus deterministic minimization 

In principle, within a stochastic approach, the exact 
minimum is never reached as the forces fi are known 
only statistically with some error bar A/j. We have found 
that the method becomes efficient when all the forces are 
non zero only within few tenths of standard deviations. 
Then for small enough constant At, and large enough bin 
compatible with the computer resources, the stochastic 
minimization, obtained with a statistical evaluation of s 
continues in a stable manner, and all the variational pa- 
rameters fluctuate after several iterations around a mean 
value. After averaging these variational parameters, the 



corresponding mean values represent very good estimates 
satisfying the minimum energy condition. This can be 
verified by performing an independent Monte Carlo sim- 
ulation much longer than the bin used for a single itera- 
tion of the stochastic minimization, and then by explic- 
itly checking that all the forces fi are zero within the 
statistical accuracy. An example is given in Fig |3 and 
discussed in the previous subsection. 

On the other hand, whenever few variational parame- 
ters are clearly out of minimum, with \fi/Afi\ > a cut , 
with a cu t — 10, we have found a faster convergence 
with a much larger At, by applying the minimization 
scheme only for those selected parameters such that 
1/,/A/il > a cut , until |/i/A/j| < a cut . After this initial- 
ization it is then convenient to proceed with the global 
minimization with all parameters changed at each itera- 
tion, in order to explore the variational space in a much 
more effective way. 



D. Different energy scales 

The SR method performs generally very well, when- 
ever there is only one energy scale in the variational wave 
function. However if there are several energy scales in the 
problem, some of the variational parameters, e.g. the 
ones defining the low energy valence orbitals, converge 
very slowly with respect to the others, and the number 
of iterations required for the equilibration becomes ex- 
ceedingly large, considering also that the time step At 
necessary for a stable convergence depends on the high 
energy orbitals, whose dynamics cannot be accelerated 
beyond a certain threshold. It is easy to understand that 
SR technique not necessarily becomes inefficient for ex- 
tensive systems with large number of atoms. Indeed sup- 
pose that we have N atoms very far apart so that we 
can neglect the interaction between electrons belonging 
to different atoms, than it is easy to see that the stochas- 
tic matrix Eq. 1141 factorizes in N smaller matrices and 
the At necessary for the convergence is equal to the sin- 
gle atom case, simply because the variational parameters 
of each single atom can evolve independently form each 
other. This is due to the size consistency of our trial func- 
tion that can be factorized as a product of N single atom 
trial functions in that limit. Anyway for system with a 
too large energy spread a way to overcome this difficulty 
was presented recently in Ref. Unfortunately this 
method is limited to the optimization of the variational 
parameters in a super-CI-basis, to which a Jastrow term 
is applied, that however can not be optimized together 
with the CI coefficients. 

In the present work, limited to a rather small atomic 
basis, the SR technique is efficient, and general enough 
to perform the simultaneous optimization of the Jastrow 
and the determinantal part of the wave function, a very 
important feature that allows to capture the most non 
trivial correlations contained in our variational ansatz. 
Moreover, SR has been extended to deal with the struc- 
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tural optimization of a chemical system, which is another 
appealing characteristic of this method. The results pre- 
sented in the next section show that in some non trivial 
cases the chemical accuracy can be reached also within 
this framework. 

However we feel that an improvement along the line 
described in Ref. \lQ will be useful for realistic electronic 
simulations of complex systems with many atoms, or 
when a very high precision is required at the variational 
level and consequently a wide spread of energy scales has 
to be included in the atomic basis. We believe that the 
difficulty to work with a large basis set will be possibly 
alleviated by using pscudopotcntials that allows to avoid 
the high energy components of the core electrons. How- 
ever more work is necessary to clarify the efficiency of 
the SR minimization scheme described here. 



IV. APPLICATION OF THE JAGP TO 
MOLECULES 

In this work we study total, correlation, and atom- 
ization energies, accompanied with the determination of 
the ground state optimal structure for a restricted en- 
semble of molecules. For each of them we performed a 
full all-electron SR geometry optimization, starting from 
the experimental molecular structure. After the energy 
minimization, we carried out all-electron VMC and DMC 
simulations at the optimal geometry within the so called 
" fixed node approximation" . The basis used here is a 
double zcta Slater set of atomic orbitals (STO-DZ) for 
the AGP part, while for parameterizing the 3-body Jas- 
trow geminal we used a double zeta Gaussian atomic set 
(GTO-DZ). In this way both the antisymmetric and the 
bosonic part are well described, preserving the right ex- 
ponential behavior for the former and the strong local- 
ization properties for the latter. Sometimes, in order 
to improve the quality of the variational wave function 
we employed a mixed Gaussian and Slater basis set in 
the Jastrow part, that allows to avoid a too strong de- 
pendency in the variational parameters in a simple way. 
However, both in the AGP and in the Jastrow sector we 
never used a large basis set, in order to keep the wave 
function as simple as possible. The accuracy of our wave 
function can be obviously improved by an extension of 
the one particle basis set. but, as discussed in the pre- 
vious section, this is rather difficult for a stochastic min- 
imization of the energy. Nevertheless, for most of the 
molecules studied with a simple JAGP wave function, a 
DMC calculation is able to reach the chemical accuracy 
in the binding energies and the SR optimization yields 
very precise geometries already at the VMC level. 

In the first part of this section some results will be 
presented for a small set of widely studied molecules and 
belonging to the Gl database. In the second subsection 
we will treat the benzene and its radical cation CqH^ , 
by taking into account its distortion due to the Jahn- 
Tcllcr effect, that is well reproduced by our SR geometry 



optimization. 



A. Small diatomic molecules, methane, and water 

Except from Bei and C2, all the molecules presented 
here belong to the standard Gl reference set: all their 
properties are well known and well reproduced by stan- 
dard quantum chemistry methods, therefore they con- 
stitute a good case for testing new approaches and new 
wave functions. 

The Li dimer is one of the easiest molecules to be 
studied after the H2, which is exact for any Diffusion 
Monte Carlo (FN DMC) calculation with a trial wave 
function that preserves the nodeless structure. L12 is 
less trivial due to the presence of core electrons that are 
only partially involved in the chemical bond and to the 
2s — 2p near degeneracy for the valence electrons. There- 
fore many authors have done benchmark calculation on 
this molecule to check the accuracy of the method or to 
determine the variance of the inter-nuclear force calcu- 
lated within a QMC framework. In this work we start 
from L12 to move toward a structural analysis of more 
complex compounds, thus showing that our QMC ap- 
proach is able to handle relevant chemical problems. In 
the case of Li2, a 3s lp STO — DZ AGP basis and a 
Is lp GTO—DZ Jastrow basis turns out to be enough for 
the chemical accuracy (see the Appendix [U] for a detailed 
description of the trial wave function). More than 99% of 
the correlation energy is recovered by a DMC simulation 
(Table P), and the atomization energy is exact within few 
thousandth of eV (0.02 kcal mo/" 1 ) (Table UH . 

Similar accuracy have been previously reached within 
a DMC approach p9j. only by using a multi-reference 
CI like wave function, that before our work, was the 
usual way to improve the electronic nodal structure. As 
stressed before, the JAGP wave function includes many 
resonating configurations through the geminal expansion, 
beyond the Is 2s HF ground state. The bond length has 
been calculated at the variational level through the fully 
optimized JAGP wave function: the resulting equilib- 
rium geometry turns out to be highly accurate (Table 
llllf) . with a discrepancy of only O.OOlcio from the ex- 
act result. For this molecule it is worth comparing our 
work with the one by Assaraf and Caffarel [3(|. Their 
zero-variance zero-bias principle has been proved to be 
effective in reducing the fluctuations related to the inter- 
nuclear force; however they found that only the inclusion 
of the space warp transformation drastically lowers the 
force statistical error, which magnitude becomes equal or 
even lower than the energy statistical error, thus allow- 
ing a feasible molecular geometry optimization. Actu- 
ally, our way of computing forces (see Eq. I27|l provides 
slightly larger variances, without explicitly invoking the 
zero- variance zero-bias principle. 

The very good bond length, we obtained, is proba- 
bly due to two main ingredients of our calculations: we 
have carried out a stable energy optimization that is of- 
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ten more effective than the variance one, as shown by 
different authors |3l| . and we have used very accurate 
trial function as it is apparent from the good variational 
energy. 

Indeed within our scheme we obtain good results with- 
out exploiting the computationally much more demand- 
ing DMC, thus highlighting the importance of the SR 
minimization described in Subsection IIII Bl 

Let us now consider larger molecules. Both Ci and O2 
are poorly described by a single Slater determinant, since 
the presence of the nondynamic correlation is strong. In- 
stead with a single geminal JAGP wave function, includ- 
ing implicitly many Slater-determinants|lQj . it is possi- 
ble to obtain a quite good description of their molecular 
properties. For C2, we used a 2s lp STO — DZ basis in 
the geminal, and a 2s lp DZ Gaussian Slater mixed basis 
in the Jastrow, for 2 we employed a 3s lp STO — DZ 
in the geminal and the same Jastrow basis as before. In 
both the cases, the variational energies recover more than 
80% of the correlation energy, the DMC ones yield more 
than 90%, as shown in Tab. [I] These results are of the 
same level of accuracy as those obtained by Filippi et 
a/|29j with a multireference wave function by using the 
same Slater basis for the antisymmetric part and a differ- 
ent Jastrow factor. ^From the TablclTTIof the atomization 
energies, it is apparent that DMC considerably improves 
the binding energy with respect to the VMC values, al- 
though for these two molecules it is quite far from the 
chemical accuracy (~ 0.1 eV): for C2 the error is 0.55(3) 
cV, for O2 it amounts to 0.67(5) eV. Indeed, it is well 
known that the electronic structure of the atoms is de- 
scribed better than the corresponding molecules if the ba- 
sis set remains the same, and the nodal error is not com- 
pensated by the energy difference between the separated 
atoms and the molecule. In a benchmark DMC calcula- 
tion with pseudopotentials [32| . Grossman found an er- 
ror of 0.27 eV in the atomization energy for O2, by using 
a single determinant wave function; probably, pseudopo- 
tentials allow the error between the pseudoatoms and the 
pseudomolecule to compensate better, thus yielding more 
accurate energy differences. As a final remark on the O2 
and C2 molecules, our bond lengths are in between the 
LDA and GGA precision, and still worse than the best 
CCSD calculations, but our results may be considerably 
improved by a larger atomic basis set, that we have not 
attempted so far. 

Methane and water arc very well described by the 
JAGP wave function. Also for these molecules we re- 
cover more than 80% of correlation energy at the VMC 
level, while DMC yields more than 90%, with the same 
level of accuracy reached in previous Monte Carlo studies 
[33l l34l l35l |3(| . Here the binding energy is almost exact, 
since in this case the nodal energy error arises essentially 
from only one atom (carbon or oxygen) and therefore it 
is exactly compensated when the atomization energy is 
calculated. Also the bond lengths arc highly accurate, 
with an error lower then 0.005 clq. 

For Be2 we applied a 3s lp STO-DZ basis set for the 



AGP part and a 2s 2p DZ Gaussian Slater mixed ba- 
sis for the Jastrow factor. VMC calculations performed 
with this trial function at the experimental equilibrium 
geometry yield 90% of the total correlation energy, while 
DMC gives 97.5% of correlation, i.e. a total energy of 
-29.33341(25) H. Although this value is better than that 
obtained by Filippi et al [29J (-29.3301(2) H) with a 
smaller basis (3s atomic orbitals not included), it is not 
enough to bind the molecule, because the binding energy 
remains still positive (0.0069(37) H). Instead, once the 
molecular geometry has been relaxed, the SR optimiza- 
tion finds a bond distance of 13.5(5) ao at the VMC level; 
therefore the employed basis allows the molecule to have 
a Van der Waals like minimum, quite far from the experi- 
mental value. In order to have a reasonable description of 
the bond length and the atomization energy, one needs to 
include at least a 3s2p basis in the antisymmetric part, as 
pointed out in Ref . 1371 and indeed an atomization energy 
compatible with the experimental result (0.11(1) eV) has 
been obtained within the extended geminal model [381 ] by 
using a much larger basis set (9s,7p,4d,2f,lg) [3!|. This 
suggests that a complete basis set calculation with JAGP 
may describe also this molecule. However, as already 
mentioned in subsection IIII Dl our SR method can not 
cope with a very large basis in a feasible computational 
time. Therefore we believe that at present the accuracy 
needed to describe correctly Be2 is out of the possibilities 
of the approach. 



B. Benzene and its radical cation 

We studied the 1 A lg ground state of the benzene 
molecule by using a very simple one particle basis set: for 
the AGP, a 2slp DZ set centered on the carbon atoms and 
a Is SZ on the hydrogen, instead for the 3-body Jastrow, 
a lslp DZ-GTO set centered only on the carbon sites. 
CqHq is a peculiar molecule, since its highly symmetric 
ground state, which belongs to the Dgh point group, is 
a resonance among different many-body states, each of 
them characterized by three double bonds between car- 
bon atoms. This resonance is responsible for the stability 
of the structure and therefore for its aromatic properties. 
We started from a non resonating 2-body Jastrow wave 
function, which dimerizes the ring and breaks the full ro- 
tational symmetry leading to the Kekule configuration. 
As we expected, the inclusion of the resonance between 
the two possible Kekule states lowers the VMC energy by 
more than 2 eV. The wave function is further improved by 
adding another type of resonance, that includes also the 
Dewar contributions connecting third nearest neighbor 
carbons. As reported in Tab. IIVI the gain with respect 
to the simplest Kekule wave function amounts to 4.2 eV, 
but the main improvement arises from the further inclu- 
sion of the three body Jastrow factor, which allows to 
recover the 89% of the total atomization energy at the 
VMC level. The main effect of the three body term is to 
keep the total charge around the carbon sites to approx- 
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imately six electrons, thus penalizing the double occupa- 
tion of the p z orbitals. The same important correlation 
ingredient is present in the well known Gutzwiller wave 
function already used for polyacetylene |4jj • Within this 
scheme we have systematically included in the 3-body 
Jastrow part the same type of terms present in the AGP 
one, namely both g a b and X a b are non zero for the same 
pairs of atoms. As expected, the terms connecting next 
nearest neighbor carbon sites are much less important 
than the remaining ones because the VMC energy does 
not significantly improve (see the full resonating + 3- 
body wave function in Tab. IIV| l. A more clear behavior 
is found by carrying out DMC simulations: the interplay 
between the resonance among different structures and 
the Gutzwiller-like correlation refines more and more the 
nodal surface topology, thus lowering the DMC energy 
by significant amounts. Therefore it is crucial to insert 
into the variational wave function all these ingredients in 
order to have an adequate description of the molecule. 
For instance, in Fig. |||we report the density surface dif- 
ference between the non-resonating 3-body Jastrow wave 
function, which breaks the C§ rotational invariance, and 
the resonating Kekule structure, which preserves the cor- 
rect A\ g symmetry: the change in the electronic structure 
is significant. The best result for the binding energy is 



the latest frozen core CCSD(T) calculations are 
able to reach a precision of 0.1 eV, but only after the 
complete basis set extrapolation and the inclusion of the 
core valence effects to go beyond the pseudopotential ap- 
proximation. Without the latter corrections, the error 
is quite lar ge e ven in the CCSD approach, amounting 
to 0.65 eV |42|. In our case, such an error arises from 
the fixed node approximation, whose nodal error is not 
compensated by the difference between the atomic and 
the molecular energies, as already noticed in the previous 
subsection. 

The radical cation CqHq of the benzene molecule has 
been the subject of intense theoretical studies^, HEj . 
aimed to focus on the Jahn- Teller distorted ground state 
structure. Indeed the ionized 2 E\ g state, which is degen- 
erate, breaks the symmetry and experiences a relaxation 
from the Dq^ point group to two different states, 2 i?2g 
and 2 B% g , that belong to the lower D2h point group. In 
practice, the former is the elongated acute deformation of 
the benzene hexagon, the latter is its compressed obtuse 
distortion. We applied the SR structural optimization, 
starting from the 2 E\ g state, and the minimization cor- 
rectly yielded a deformation toward the acute structure 
for the 2 B 2g state and the obtuse for the 2 B 3g one; the 
first part of the evolution of the distances and the angles 
during those simulations is shown in Fig|SJ After this 
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FIG. 4: Surface plot of the charge density projected onto the 
molecular plane. The difference between the non-resonating 
(indicated as HF) and resonating Kekule 3-body Jastrow wave 
function densities is shown. Notice the corresponding change 
from a dimerized structure to a Ce rotational invariant density 
profile. 

obtained with the Kekule Dcwar resonating 3 body wave 
function, which recovers the 98,6% of the total atom- 
ization energy with an absolute error of 0.84(8) eV. As 
Pauling |4l| first pointed out, benzene is a genuine RVB 
system, indeed it is well described by the JAGP wave 
function. Moreover Pauling gave an estimate for the res- 
onance energy of 1.605 eV from thermochemical experi- 
ments in qualitative agreement with our results. A final 
remark about the error in the total atomization energy: 
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FIG. 5: Plot of the convergence toward the equilibrium geom- 
etry for the 2 B2 g acute and the 2 Bz g obtuse benzene cation. 
Notice that both the simulations start form the ground state 
neutral benzene geometry and relax with a change both in 
the C — C bond lengths and in the angles. The symbols are 
the same of Tab. 

equilibration, average over 200 further iterations yields 
bond distances and angles with the same accuracy as 
the all-electron BLYP/6-31G* calculations reported in 
Ref. (see Tab. E). As it appears from Tab. EU 
not only the structure but also the DMC total energy is 
in perfect agreement with the BLYP/6-31G*, and much 
better than SVWN/6-31G* that does not contain semi 
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empirical functional, for which the comparison with our 
calculation is more appropriate, being fully ab-initio. 

The difference of the VMC and DMC energies between 
the two distorted cations are the same within the error 
bars; indeed, the determination of which structure is the 
real cation ground state is a challenging problem, since 
the experimental results give a difference of only few meV 
in favor of the obtuse state and also the most refined 
quantum chemistry methods are not in agreement among 
themselves |44j- A more affordable problem is the de- 
termination of the adiabatic ionization potential (A1P), 
calculated for the 2 B^ g state, following the experimental 
hint. Recently, very precise CCSD(T) calculations have 
been performed in order to establish a benchmark theo- 
retical study for the ionization threshold of benzene ; 
the results are reported in Tab. lVffl After the correction 
of the zero point energy due to the different structure of 
the cation with respect to the neutral molecule and taken 
from a B3LYP / cc-pVTZ calculation reported in Ref . 0, 
the agreement among our DMC result, the benchmark 
calculation and the experimental value is impressive. No- 
tice that in this case there should be a perfect cancella- 
tion of nodal errors in order to obtain such an accurate 
value; however, we believe that it is not a fortuitous re- 
sult, because in this case the underlying nodal structure 
does not change much by adding or removing a single 
electron. Therefore we expect that this property holds 
for all the affinity and ionization energy calculations with 
a particularly accurate variational wave function as the 
one we have proposed here. Nevertheless DMC is needed 
to reach the chemical accuracy, since the VMC result is 
slightly off from the experimental one just by few tenths 
of eV. The AIP and the geometry determination for the 
CqHq arc encouraging to pursue this approach, with the 
aim to describe even much more interesting and challeng- 
ing chemical systems. 



V. CONCLUSION 

In this work, we have tested the JAGP wave function 
on simple molecular systems where accurate results are 
known. As shown in the previous section a large amount 
of the correlation energy is already recovered at the varia- 
tional level with a computationally very efficient and fea- 
sible method, extended in this work to the nuclear geom- 
etry optimization. Indeed, much larger systems should 
be tractable because, within the JAGP ansatz, it is suf- 
ficient to sample a single determinant whose dimension 
scales only with the number of electrons. The presence 
of the Jastrow factor implies the evaluation of multi- 
dimensional integrals that, so far can be calculated ef- 
ficiently only with the Monte Carlo method. Within this 
framework, it is difficult to reach the complete basis set 
limit, both in the Jastrow and the AGP terms, although 
some progress has been made recentlv|l9t l46j . Even if 
the dimension of the basis is limited by the difficulty to 
perform energy optimization with a very large number 



of variational parameters, we have obtained the chemical 
accuracy for most cases studied. From a general point of 
view the basis set convergence of the JAGP is expected to 
be faster than AGP considering that the electron-electron 
cusp condition is fulfilled exactly at each level of the ex- 
pansion. Nevertheless, all results presented here can be 
systematically improved with larger basis set. In partic- 
ular the Be2 bonding distance should be substantially 
corrected by a more complete basis, that we have not 
attempted so far.[39l| 

The usefulness of the JAGP wave function is already 
well known in the study of strongly correlated systems 
defined on a lattice. For instance in the widely studied 
Hubbard model, as well as in any model with electronic 
repulsion, it is not possible to obtain a superconducting 
ground state at the mean-field Hartree-Fock level. On the 
contrary, as soon as a correlated Jastrow term is applied 
to the BCS wave function (equivalent to the AGP wave 
function in momentum spacepTfl'). the stabilization of a 
d-wave superconducting order parameter is possible, and 
is expected to be a realistic property of the model jlgj. 
More interestingly the presence of the Jastrow factor can 
qualitatively change the wave function especially at one 
electron per site filling, by converting a BCS supercon- 
ductor to a Mott insulator with a finite charge gap^ij- 

The same effect is clearly seen for the gedanken experi- 
ment of a dilute gas of H 2 molecules, a clarifying test ex- 
ample already used in the introduction. The AGP wave 
function is essentially exact for a single molecule (at least 
with the complete basis set), but its obvious size consis- 
tent extension to the gas would lead to the unphysical 
result of superconductivity because the charge around 
each molecule would be free to fluctuate within the cho- 
sen set of geminal orbitals. Only the presence of the Jas- 
trow term added to this wave function, allows the local 
conservation of the charge around the molecule, by for- 
bidding unphysical H2 dimers with more than two elec- 
trons. Once the charge is locally conserved, the phase of 
the BCS- AGP wave function cannot have a definite value 
and phase coherence is correctly forbidden by the Jas- 
trow factor. In the present work, the interplay between 
the Jastrow and the geminal part has been shown to be 
very effective in all cases studied and particularly in the 
non trivial case of the benzene molecule, where we have 
shown systematically the various approximations. Only 
when both the Jastrow and the AGP terms are accurately 
optimized together, the AGP nodal structure of the wave 
function is considerably improved. For the above reasons 
and the size consistency of the JAGP we expect that this 
wave function should be generally accurate also in com- 
plex systems made by many molecules. The local con- 
servation of the charge around each molecule is taken 
into account by the Jastrow factor, whereas the quality 
of each molecule is described also by the AGP geminal 
part exactly as in the H2 gas example. 

In the near future it is very appealing and promising to 
extend the JAGP study to the DNA nitrogenous bases, 
whose geometrical structure is very similar to the benzene 
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ring. In particular we plan to accurately evaluate the en- 
ergetics (reduction potential, ionization energies, electron 
affinity, etc.) of DNA bases and base pairs, quantities of 
great importance to characterize excess electron and hole 
transfer which are involved in radiation damage as well 
as in the development of novel DNA technologies. 



APPENDIX A: STABILIZATION OF THE SR 
TECHNIQUE 

Whenever the number of variational parameters in- 
creases, it often happens that the stochastic (N + 1) x 
(N + 1) matrix 



determinant S of the corresponding Nx N overlap matrix 



Sk,k< — 



(*\OkOk>\*) 

(*l*> 



(Al) 



becomes singular, i.e. the condition number, defined as 
the ratio a = Aat/Ai between its maximum A at and min- 
imum eigenvalue Ai, is too large. In that case the inver- 
sion of this matrix generates clear numerical instabilities 
which are difficult to control especially within a statisti- 
cal method. Here Ok ~ dh ^ x ^ arc the operators corre- 
sponding to the variational parameters ak appearing in 
the wave function \& for k = 1, • • • N, whereas for k = 
the operator Oq represents the identity one. 

The first successful proposal to control this instabil- 
ity was to remove from the inversion problem (|12f) . re- 
quired for the minimization, those directions in the vari- 
ational parameter space corresponding to exceedingly 
small eigenvalues A^. 

In this appendix we describe a better method. As a 
first step, we show that the reason of the large condition 
number a is due to the existence of "redundant" varia- 
tional parameters that do not make changes to the wave 
function within a prescribed tolerance e. Indeed in prac- 
tical calculations, we are interested in the minimization of 
the wave function within a reasonable accuracy. The tol- 
erance e may represent therefore the distance between the 
exact normalized variational wave function which mini- 
mizes the energy expectation value and the approximate 
acceptable one, for which we no longer iterate the min- 
imization scheme. For instance e = 1/1000 is by far 
acceptable for chemical and physical interest. A stable 
algorithm is then obtained by simply removing the pa- 
rameters that do not change the wave function by less 
than e from the minimization. An efficient scheme to 
remove the "redundant parameters" is also given. 

Let us consider the N normalized states orthogonal to 
^P, but not orthogonal among each other: 



(Ofc-«fc,o)l*) 



(A2) 



where s^g is defined in Ea. lAll These normalized vectors 
define N directions in the N— dimensional variational pa- 
rameter manifold, which are independent as long as the 



Sk.k' 



(A3) 



is non zero. The number S is clearly positive and it as- 
sumes its maximum value 1 whenever all the directions 
ei are mutually orthogonal. On the other hand, let us 
suppose that there exists an eigenvalue A of s smaller 
than the square of the desired tolerance e 2 , then the cor- 
responding eigenvector \v >= J2i a i\ e i) ^ s sucn that: 



A 



(A4) 



where the latter equation holds due to the normalization 
condition J2i a i — 1- We arrive therefore to the conclu- 
sion that it is possible to define a vector v with almost 
vanishing norm \v\ = \f\ < e as a linear combination of 
ei, with at least some non zero coefficient. This implies 
that the N directions ek are linearly dependent within 
a tolerance e and one can safely remove at least one pa- 
rameter from the calculation. 

In general whenever there are p vectors v.i that are 
below the tolerance e the optimal choice to stabilize 
the minimization procedure is to remove p rows and p 
columns from the matrix l|A3[) . in such a way that the 
corresponding determinant of the (N — p) x (N — p) over- 
lap matrix is maximum. 

^From practical purposes it is enough to consider an 
iterative scheme to find a large minor, but not necessarily 
the maximum one. This method is based on the inverse of 
s. At each step we remove the i — th row and column from 
s for which s~l is maximum. We stop to remove rows and 
columns after p inversions. In this approach we exploit 
the fact that, by a consequence of the Laplace theorem 
on determinants, s^], is the ratio between the described 
minor without the k — th row and column and the de- 
terminant of the full s matrix. Since within a stochastic 
method it is certainly not possible to work with a machine 
precision tolerance, setting e = 0.001 guarantees a stable 
algorithm, without affecting the accuracy of the calcula- 
tion. The advantage of this scheme, compared with the 
previous one|l7j. is that the less relevant parameters can 
be easily identified after few iterations and do not change 
further in the process of minimization. 



APPENDIX B: SIZE CONSISTENCY OF THE 
3-BODY JASTROW FACTOR 

In order to prove the size consistency property of the 
three body Jastrow factor described in Sec. Ill CI let us 
take into account a system composed by two well sepa- 
rated subsystems A and B, which are distinguishable and 
whose dimensions are much smaller than the distance be- 
tween themselves; in general they may contain more then 
one atom. In this case the Jastrow function J3 l|10|) can 
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be written as J3 = e u with: 



rearranged in the following form: 



i,jeA iJeB 

ieA j£B 

where we have explicitly considered the sum over dif- 
ferent subsystems. As usual, the two particle function 
4 > ( r ii r j) is expanded over a single particle basis ip, cen- 
tered on each nucleus of the system: 

<l>{r i ,r j )=Y,X n ' n ^ m {riW{r j ). (B2) 

The indices n and m refer not only to the basis elements 
but also to the nuclei which the orbitals are centered on. 

The self consistency problem arises from the last term 
in Ea. lBll i.e. when the electron r$ belongs to A and rj 
to B. If the Jastrow is size consistent, whenever A and B 
are far apart from each other this term must vanish or at 
most generate a one body term that is in turn size con- 
sistent, as we are going to show in the following. In the 
limit of large separation all the \ m n off diagonal terms 
connecting any basis element of A to any basis element 
of B must vanish. The second requirement is a suffi- 
ciently fast decay of the basis set orbitals ip(r) whenever 
r — > 00, except at most for a constant term C n which may 
be present in the single particle orbitals, and is useful to 
improve the variational energy. 

For the sake of generality, suppose that the system A 
contains Ma nuclei and Na electrons. The first require- 
ment implies that: 



(n,rj) = A m '"^ m (r i )^"(r i ) (B3) 

+ ^ m ' n r\n)r(r 3 ), 



instead the second allows to write the following expres- 
sion for the mixed term in Eq. IB1I 

i£Aj£B n£A m£B 

(B4) 

where the factors P n are one body terms defined as: 

Notice that if all the orbitals decay to zero, the size 
consistency is immediately recovered, since the sum in 
Eg. IB4I vanishes. Analogously to the derivation we have 
done to extract the one body contribution from the mixed 
term, the other two terms on the RHS of Eg. IB II can be 



\ J2 Mn,^) = (N A -l)J2CnPn (B6) 

i,j€A n&A 
i ¥=3 

+ two body terms, 
and the sum in Eq. IB II can be rewritten as: 

U = (N-1)J2 CnPn + (N-1)J2 CnPn (B7) 
neA neB 

+ two body size consistent terms. 

Therefore the size consistency implies that the scaling of 
the C n with the total number of particle N is: 



C n 



as mentioned in Sec. Ill CI 



N - 1 



(B8) 



APPENDIX C: AN EXAMPLE CASE: 
WAVE FUNCTION FOR Li 2 



JAGP 



We briefly describe the application of the JAGP to 
the L%2 molecule. This example shows the beauty of our 
approach that allows to describe the chemical bond as 
resonance of many pairing functions whose importance 
is weighted by the A coefficients. In the expansion of the 
geminal function for the determinant in Eq. we used 
six orbitals for each atom: 



Hs,(P2s 
fop 
03s 



= C U 



■ pe 



C 3s r 2 e-^ r . 



(CI) 
(C2) 
(C3) 



The parameters p in Is, 2s orbitals are fixed by the single 
atomic cusp conditions, and C\ s ,C2p,C^ s are the nor- 
malization constants. These orbitals are connected by 
different A^ h „, which obey the constraints given by the 
symmetry of the system, and are reported in table IVIlTl 
Since the trial function does not need to be normalized 
we set X\'si s equal to 1. The total number of non zero 
A is 58, but the constraints allow to reduce them to 18 
variational parameters. 

In the Jastrow part we used a two body term that is a 
slightly modified version of the Eq. [5J In fact due to the 
simple symmetry of the system is possible to build a Jas- 
trow more suitable for this diatomic molecule, namely: 



u(r, z) 



2(1+ v/aOr 2 + y 2 ) + 6z 2 ) ' 



(C4) 



which distinguishes the different components of the two 
electrons distance. We found that this two body Jas- 
trow factor is particularly useful for Li 2 , which is much 
more elongated than the other molecules studied here, 
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for which the usual form in Eq. has been employed. 
The optimal parameters obtained for the Jastrow are 
a = 0.8796 , b — 0.7600. In the expansion of the pairing 
function for the three body Jastrow term (see Eq. I10|) 
we used the following orbitals: 

0sG = e~'^+p, (C5) 
<t> vG = r(e~^ 2 + P e- Z ^ . (C6) 

The A matrix that connects these orbitals is given in table 
IIXI this matrix fulfills the same symmetry constraints 
as in the case of the paring determinant. In this case 
the total number of non zero A is 24 and the symmetry 



reduces the variational freedom to only 8 parameters. 
The single particle orbitals are reported in table Ixl and 
include other 15 parameters. 
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TABLE I: Total energies in variational (Evmc) and diffusion (Edmc) Monte Carlo calculations; the percentages of correlation 
energy recovered in VMC (EY MC (%)) and DMC (E° MC (%)) have been evaluated using the "exact" (E ) and Hartree-Fock 
(Ehf) energies from the references reported in the table. Here "exact" means the ground state energy of the non relativistic 
infinite nuclear mass hamiltonian. The energies are in Hartree. 





E 


Ehf 


Evmc 


eY mc (%) 


Edmc 


E? MC (%) 


Li 


-7.47806° 


-7.432727° 


-7.47721(11) 


98.12(24) 


-7.47791(12) 


99.67(27) 




-14.9954 c 


-14.87152 c 


-14.99002(12) 


95.7(1) 


-14.99472(17) 


99.45(14) 


Be 


-14.66736 a 


-14.573023° 


-14.66328(19) 


95.67(20) 


-14.66705(12) 


99.67(13) 


Be 2 


-29.33854(5) c 


-29.13242 c 


-29.3179(5) 


89.99(24) 


-29.33341(25) 


97.51(12) 


O 


-75.0673° 


-74.809398° 


-75.0237(5) 


83.09(19) 


-75.0522(3) 


94.14(11) 


H 2 


-76.438(3)° 


-76.068(1)° 


-76.3803(4) 


84.40(10) 


-76.4175(4) 


94.46(10) 


O2 


-150.3268 c 


-149.6659 c 


-150.1992(5) 


80.69(7) 


-150.272(2) 


91.7(3) 


C 


-37.8450° 


-37.688619° 


-37.81303(17) 


79.55(11) 


-37.8350(6) 


93.6(4) 


c 2 


-75.923(5) c 


-75.40620 c 


-75.8273(4) 


81.48(8) 


-75.8826(7) 


92.18(14) 


CH4 


-40.515 d 


-40.219 d 


-40.4627(3) 


82.33(10) 


-40.5041(8) 


96.3(3) 


C 6 H 6 


-232.247(4) c 


-230.82(2) / 


-231.8084(15) 


69.25(10) 


-232.156(3) 


93.60(21) 



"Exa ct a nd HF energies from Ref . l5Ct 

kRef.py 

c Ref. 2C 

d Ref.|3S 

e Esti mat ed "exact" energy from Ref. |4 1 
^Ref.H3 



TABLE II: Binding energies in eV obtained by variational (Avmc) and diffusion (Admc) Monte Carlo calculations; Ao is 
the "exact" result for the non-relativistic infinite nuclear mass hamiltonian. Also the percentages (Avmc(%) and Ada/c(%)) 
of the total binding energies are reported. 



Ao Avmc Avmc{%) Admc A_da/c(%) 

Li 2 -1.069 -0.967(3) 90.4(3) -1.058(5) 99.0(5) 

O2 -5.230 -4.13(4) 78.9(8) -4.56(5) 87.1(9) 

H 2 -10.087 -9.704(24) 96.2(1.0) -9.940(19) 98.5(9) 

C* 2 -6.340 -5.476(11) 86.37(17) -5.79(2) 91.3(3) 

CH 4 -18.232 -17.678(9) 96.96(5) -18.21(4) 99.86(22) 

C 6 H 6 -59.25 -52.53(4) 88.67(7) -58.41(8) 98.60(13) 



TABLE III: Bond lengths (R) in atomic units; the subscript refers to the "exact" results. For the water molecule R is the 
distance between O and H and is the angle HOH (in deg), for CH4 R is the distance between C and H and 9 is the HCH 
angle. 





Ro 


R 


0o 


9 


Li 2 


5.051 


5.0516(2) 






O2 


2.282 


2.3425(18) 






C2 


2.348 


2.366(2) 






H 2 


1.809 


1.8071(23) 


104.52 


104.74(17) 


CH4, 


2.041 


2.049(1) 


109.47 


109.55(6) 




nCC 


R ua 


nCH 


R UH 


Cells 


2.640 


2.662(4) 


2.028 


1.992(2) 
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TABLE IV: Binding energies in eV obtained by variational (Avaic) and diffusion (Admc) Monte Carlo calculations with 
different trial wave functions for benzene. In order to calculate the binding energies yielded by the 2-body Jastrow we used 
the atomic energies reported in Ref. llOt The percentages (Avmc{%) and Admc(%)) of the total binding energies are also 
reported. 



Kekule + 2body 
resonating Kekule + 2body 
resonating Dewar Kekule + 2body 
Kekule + 3body 
resonating Kekule + 3body 
resonating Dewar Kekule + 3body 
full resonating + 3body 



-(%) 



Ai 



-30.57(5) 
-32.78(5) 
-34.75(5) 
-49.20(4) 
-51.33(4) 
-52.53(4) 
-52.65(4) 



51.60(8) 
55.33(8) 
58.66(8) 
83.05(7) 
86.65(7) 
88.67(7) 
88.87(7) 



-56.84(11) 

-55.54(10) 

-57.25(9) 

-58.41(8) 

-58.30(8) 



95.95(18) 
93.75(17) 
96.64(15) 
98.60(13) 
98.40(13) 



TABLE V: Bond lengths (r) for the two lowest 2 B 2g and 2 B 3g states of the benzene radical cation. The angles a are expressed 
in degrees, the lengths in do. The carbon sites are numerated from 1 to 6. 





2 B 2 g 

acute 


2 B 3g 
obtuse 


Computational method 


r(Ci - C* 2 ) 


2.616 


2.694 


B3LYP/cc-pVTZ " 




2.649 


2.725 


BLYP/6-31G* b 




2.659(1) 


2.733(4) 


SR-VMC c 


r(C 2 - C 3 ) 


2.735 


2.579 


B3LYP/cc-pVTZ 




2.766 


2.615 


BLYP/6-31G* b 




2.764(2) 


2.628(4) 


SR-VMC c 


ol(C§CiC2) 


118.4 


121.6 


B3LYP/cc-pVTZ ° 




118.5 


121.5 


BLYP/6-31G* b 




118.95(6) 


121.29(17) 


SR-VMC c 



i, Ref.E3 

c This work 



TABLE VI: Total energies for the 2 B 2g and 2 B 3g states of the benzene radical cation after the geometry relaxation. A 
comparison with a BLYP/6-31G* and SVWN/6-31G* all-electron calculation (Ref.0) is reported. 

VMC PMC BLYP/6-31G* SVWN/6-31G*" 

2 B 2g -231.4834(15) -231.816(3) -231.815495 -230.547931 

2 B :ig -231.4826(14) -231.812(3) -231.815538 -230.547751 



TABLE VII: Adiabatic ionization potential of the benzene molecule; our estimate is done for the 2 B 3g relaxed geometries of 
the benzene radical cation, with an inclusion of the zero point motion correction between the 2 B 3g state and the 1 A\ g neutral 
molecule ground state, calculated in R.ef. l45l at the B3LYP/6-31G* level. 





VMC ° 


DMC a 


CCSD(T)/cc-pVooZ b 


experiment c 


AIP 


8.86(6) 


9.36(8) 


9.29(4) 




AZPE ad 


-0.074 


-0.074 


-0.074 




best estimate 


8.79(6) 


9.29(8) 


9.22(4) 


9.2437(8) 



"Thi s wo rk 
c Ref.g3 
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TABLE VIII: Matrix of the A coefficients of the geminal function expansion in the pairing determinant for the L12 molecule. 
The matrix is symmetric to have a spin singlet, therefore we show only the upper part of it. 





ls a 


2s a 


2pz a 




2px a 


2py a 


3s a 




2s h 


2pz h 


2px h 


2py h 


3s b 




1 


-2.162 10~ A 


-6.838 10 


-3 








-2.877 10~ J 




















2s a 




6.22 10~ 4 


-1.601 10 


-4 








2.031 10 -3 





6.79 10~ 3 


-4.67 10 -4 








1.790 10 ~~ 3 


2p*a 






-4.25 10" 


-4 








-1.132 10~ 3 





-4.67 10 ~ 4 


3.211 10 ~ 4 








7.67 10~ 4 












-1.351 10~ 3 

















-1.173 10~ 3 




















-1.351 10~ 3 

















-1.173 10 -3 





3s a 














-1.541 10~ 3 





1.790 10 -3 


7.67 10 -4 








-8.91 10~ 4 


ls b 
















1 


-2.162 10 -3 


-6.838 10 -3 








-2.877 10~ 3 


2* 6 


















6.22 10 -4 


-1.601 10 -4 








1.790 10 ~~ 3 


2pz b 




















-4.25 10~ 4 








-1.132 10~ 3 


2px b 






















-1.351 10 -3 








2 PVb 
























-1.351 10 -3 






TABLE IX: Matrix of the A coefficients of the pairing function expansion in the three body Jastrow for the Li2 molecule. As 
in the previous table only the upper part is reported . 



sG a pGxa pGya pGz a sGt pGxb pGyb pGzb 



-0.2427 -2.713 10~ 4 -5.136 10 -4 -1.202 10 _a 

-0.1772 -7.997 10" 3 

-0.1772 -7.997 10~ 3 

1.202 10" 5 -8.749 10" 3 

-0.2427 -2.713 10" 4 

-0.1772 

-0.1772 

1.027 10 -2 



TABLE X: Orbital basis set parameters used for the L12 
molecule. Since the molecule is homonuclear the parameters 
of the atom b are the same as the atom a. 





Z\ 


Z2 


P 




2.4485 


4.2891 


0.4278 




0.5421 


1.4143 


-1.5500 


4>2px a 


0.6880 






<t>2pya 


0.6880 






fopza 


1.0528 








0.6386 






<^sG a 


1.4356 




-0.2044 


4>pGx a 


0.7969 


4.4217 


-1.2689 


4>pG Va 


0.7969 


4.4217 


-1.2689 


$pGz a 


8.980 10~ 3 


-0.1924 


0.3229 



